###########################################################
### Disorganized Political Violence:                    ###
### A Demonstration Case of Temperature and Insurgency  ###
### Replication Code:                                   ###
### Temperature Response Function and Deviation Results ### 
### Median Temperatures Robustness Check                ###
### Andrew Shaver, Alex Bollfrass                       ###
### Date: 07/04/2022                                    ###
###########################################################

# Note: these results were generated using:
# 1) MacBook Pro (16-inch, 2021); Chip: Apple M1 Max; Memory: 64 GB; running
# macOS Monterey.
# 2) R 4.1.2 GUI 1.77 High Sierra build (8007)
# 3) RStudio 2022.07.1+554 "Spotted Wakerobin" Release (7872775ebddc40635780ca1ed238934c3345c5de, 2022-07-22) for macOS
# Mozilla/5.0 (Macintosh; Intel Mac OS X 12_1_0) AppleWebKit/537.36 (KHTML, like Gecko) QtWebEngine/5.12.10 Chrome/69.0.3497.128 Safari/537.36
# (Please see below for specific package versions.)

# Load required packages:
library(MASS)
library(lubridate)
library(DataCombine)
library(broom)
library(data.table)
library(radiant.model)
library(lfe)
library(ggplot2)
library(ggpubr)
library(stargazer)
library(lfe)
library(arm)
library(locfit)

packageVersion("MASS")
# ‘7.3.55’
packageVersion("lubridate")
# ‘1.8.0’
packageVersion("DataCombine")
# ‘0.2.21’
packageVersion("broom")
# ‘0.7.12’
packageVersion("data.table")
# ‘1.14.2’
packageVersion("radiant.model")
# ‘1.4.4’
packageVersion("lfe")
# ‘2.8.7.1’
packageVersion("ggplot2")
# ‘3.3.6’
packageVersion("ggpubr")
# ‘0.4.0’
packageVersion("stargazer")
# ‘5.2.3’
packageVersion("lfe")
# ‘2.8.7.1’
packageVersion("arm")
# ‘1.12.2’
packageVersion("locfit")
# ‘1.5.9.5’

# Load data:
# For primary temperature-violence analysis:
# Daily analysis:
iq_sig1 <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/iraq_baghdad_basrah_processed_data.rds")
# For median analysis:
iq_median_temp <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/iq_daily_median.rds")

# Match median temperature data to the Iraq panel:
colnames(iq_median_temp)[2] <- "governorate"
iq_sig1 <- plyr::join(iq_sig1, iq_median_temp, by = c("date", "governorate"))

# Next, construct temperature deviations measure:
# Write function:
deviation <- function(x){
  admin_list <- list()
  library(locfit)
  library(data.table)
  for(i in 1:length(unique(panel[, x]))){
    temp <- panel[panel[, x] %in% unique(panel[, x])[i], ]
    temp$nday <- 1:nrow(temp)
    t <- locfit(median ~ lp(nday, nn = .05), data=temp)
    temp$temp_hat <- predict(t, newdata = temp$nday)
    temp$temp_dev <- temp$median - temp$temp_hat
    admin_list[[i]] <- temp
    rm(temp)
  }
  admin_list <- rbindlist(admin_list)
  return(admin_list)
}
panel <- iq_sig1
iq_sig1 <- deviation("governorate")
iq_sig1 <- as.data.frame(iq_sig1)

deviation <- function(x){
  admin_list <- list()
  library(locfit)
  library(data.table)
  for(i in 1:length(unique(panel[, x]))){
    temp <- panel[panel[, x] %in% unique(panel[, x])[i], ]
    temp$nday <- 1:nrow(temp)
    t <- locfit(median ~ lp(nday, nn = .05), data=temp)
    temp$median_hat <- predict(t, newdata = temp$nday)
    temp$median_dev <- temp$median - temp$median_hat
    admin_list[[i]] <- temp
    rm(temp)
  }
  admin_list <- rbindlist(admin_list)
  return(admin_list)
}
panel <- iq_sig1
iq_sig1 <- deviation("governorate")

iq_sig1 <- as.data.frame(iq_sig1)

# Next, generate temperature deviations for each context:
# Set temperature range for all:
temperature_range <- 65:95

###
### Deviation results:
###

# Baghdad and Basrah:
# Unconstrained Violence:
# OLS model results:

temperature_deviation_iq_ols <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$median > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(log(unconstrained+1) ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                            constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                            unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            week + governorate | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(unconstrained+1) ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                             unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + governorate | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(unconstrained+1) ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained |
                                 week + governorate | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(unconstrained+1) ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained |
                                  month + governorate | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_ols <- temperature_deviation_iq_ols(iq_sig1, 0.95)
iq_90_ols <- temperature_deviation_iq_ols(iq_sig1, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
iq_95_ols$change <- (sd(iq_sig1$median_dev, na.rm = T) * iq_95_ols$estimate) / mean(log(iq_sig1$unconstrained+1))

# For plotting Baghdad-Basrah results:
temp_ranges <- unique(iq_95_ols$bin)[c(TRUE, FALSE)]
for(i in seq(1,(length(temp_ranges)*2-2), by=2)){
  temp_ranges <- append(temp_ranges, " ", after=i)
}

# Plot:
p_iq_ols <- ggplot() + 
  
  geom_line(data = iq_95_ols[iq_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = iq_95_ols[iq_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols[iq_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols[iq_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_ols[iq_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols[iq_95_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols[iq_90_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols[iq_90_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols[iq_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols[iq_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols[iq_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols[iq_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols[iq_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols[iq_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols[iq_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols[iq_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols[iq_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = iq_95_ols[iq_95_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = iq_90_ols[iq_90_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = iq_90_ols[iq_90_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "blue4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.02, 0.0275) + 
  annotate(geom = "text", x = 65:95, y = -0.0165,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

# Negative binomial model results:

# Calculate magnitude in terms of percentage deviation from the mean:
new_data_iq <- data.frame(median_dev = c(0, sd(iq_sig1$median_dev)),
                          threshold = 0,
                          
                          prcp = mean(iq_sig1$prcp, na.rm = T), 
                          wdsp = mean(iq_sig1$wdsp, na.rm = T), 
                          dewp = mean(iq_sig1$dewp, na.rm = T), 
                          visib = mean(iq_sig1$visib, na.rm = T), 
                          mxspd = mean(iq_sig1$mxspd, na.rm = T), 
                          daylight = mean(iq_sig1$daylight, na.rm = T), 
                          hoursofpower = mean(iq_sig1$hoursofpower, na.rm = T), 
                          
                          idf = mean(iq_sig1$idf), 
                          ied = mean(iq_sig1$ied),  
                          constrained = mean(iq_sig1$constrained),  
                          unconstrained = mean(iq_sig1$unconstrained),  
                          
                          temp_1 = mean(iq_sig1$temp_1, na.rm = T),
                          temp_2 = mean(iq_sig1$temp_2, na.rm = T),
                          temp_3 = mean(iq_sig1$temp_3, na.rm = T),
                          temp_4 = mean(iq_sig1$temp_4, na.rm = T),
                          temp_5 = mean(iq_sig1$temp_5, na.rm = T),
                          temp_6 = mean(iq_sig1$temp_6, na.rm = T),
                          temp_7 = mean(iq_sig1$temp_7, na.rm = T),
                          
                          unconstrained_1 = mean(iq_sig1$unconstrained_1, na.rm = T),
                          unconstrained_2 = mean(iq_sig1$unconstrained_2, na.rm = T),
                          unconstrained_3 = mean(iq_sig1$unconstrained_3, na.rm = T),
                          unconstrained_4 = mean(iq_sig1$unconstrained_4, na.rm = T),
                          unconstrained_5 = mean(iq_sig1$unconstrained_5, na.rm = T),
                          unconstrained_6 = mean(iq_sig1$unconstrained_6, na.rm = T),
                          unconstrained_7 = mean(iq_sig1$unconstrained_7, na.rm = T),
                          
                          constrained_1 = mean(iq_sig1$constrained_1, na.rm = T),
                          constrained_2 = mean(iq_sig1$constrained_2, na.rm = T),
                          constrained_3 = mean(iq_sig1$constrained_3, na.rm = T),
                          constrained_4 = mean(iq_sig1$constrained_4, na.rm = T),
                          constrained_5 = mean(iq_sig1$constrained_5, na.rm = T),
                          constrained_6 = mean(iq_sig1$constrained_6, na.rm = T),
                          constrained_7 = mean(iq_sig1$constrained_7, na.rm = T),
                          
                          idf_1 = mean(iq_sig1$idf_1, na.rm = T),
                          idf_2 = mean(iq_sig1$idf_2, na.rm = T),
                          idf_3 = mean(iq_sig1$idf_3, na.rm = T),
                          idf_4 = mean(iq_sig1$idf_4, na.rm = T),
                          idf_5 = mean(iq_sig1$idf_5, na.rm = T),
                          idf_6 = mean(iq_sig1$idf_6, na.rm = T),
                          idf_7 = mean(iq_sig1$idf_7, na.rm = T),
                          
                          ied_1 = mean(iq_sig1$ied_1, na.rm = T),
                          ied_2 = mean(iq_sig1$ied_2, na.rm = T),
                          ied_3 = mean(iq_sig1$ied_3, na.rm = T),
                          ied_4 = mean(iq_sig1$ied_4, na.rm = T),
                          ied_5 = mean(iq_sig1$ied_5, na.rm = T),
                          ied_6 = mean(iq_sig1$ied_6, na.rm = T),
                          ied_7 = mean(iq_sig1$ied_7, na.rm = T),
                          
                          week = median(iq_sig1$week),
                          month = median(iq_sig1$month),
                          governorate = "Baghdad")

temperature_deviation_iq_nb <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$median > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week.nb <- glm(unconstrained ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                              constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                              unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                              as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint_robust(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint_robust(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_iq[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_iq[1,], type = "response")) /  mean(iq_sig1$unconstrained)
    
    reg_temp_month.nb <- glm(unconstrained ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                               temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                               unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint_robust(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint_robust(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(unconstrained ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                                   as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint_robust(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint_robust(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(unconstrained ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                                    as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint_robust(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint_robust(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_nb <- temperature_deviation_iq_nb(iq_sig1, 0.95)
iq_90_nb <- temperature_deviation_iq_nb(iq_sig1, 0.90)

# Plot:

p_iq_nb <- ggplot() + 
  
  geom_line(data = iq_95_nb[iq_95_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = iq_95_nb[iq_95_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb[iq_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb[iq_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_nb[iq_95_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb[iq_95_nb$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb[iq_90_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb[iq_90_nb$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb[iq_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb[iq_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb[iq_90_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb[iq_90_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb[iq_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb[iq_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb[iq_90_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb[iq_90_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb[iq_95_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = iq_95_nb[iq_95_nb$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = iq_90_nb[iq_90_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = iq_90_nb[iq_90_nb$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "blue4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.035, 0.05) + 
  annotate(geom = "text", x = 65:95, y = -0.03,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

# Constrained Violence:
# OLS model results:
temperature_deviation_iq_ols_con <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$median > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(log(constrained+1) ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                            constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                            unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            week + governorate | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(constrained+1) ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                             unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + governorate | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(constrained+1) ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained |
                                 week + governorate | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(constrained+1) ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained |
                                  month + governorate | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_ols_con <- temperature_deviation_iq_ols_con(iq_sig1, 0.95)
iq_90_ols_con <- temperature_deviation_iq_ols_con(iq_sig1, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
iq_95_ols_con$change <- (sd(iq_sig1$median_dev, na.rm = T) * iq_95_ols_con$estimate) / mean(log(iq_sig1$constrained+1))

# Plot:
p_iq_ols_con <- ggplot() + 
  
  geom_line(data = iq_95_ols_con[iq_95_ols_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "green4") +
  geom_line(data = iq_95_ols_con[iq_95_ols_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_con[iq_95_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_con[iq_95_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_ols_con[iq_95_ols_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_con[iq_95_ols_con$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_con[iq_90_ols_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_con[iq_90_ols_con$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_con[iq_95_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_con[iq_95_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_con[iq_90_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_con[iq_90_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_con[iq_95_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_con[iq_95_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_con[iq_90_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_con[iq_90_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_con[iq_95_ols_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "green4") + 
  geom_errorbar(data = iq_95_ols_con[iq_95_ols_con$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "green4", linetype = "dashed") +
  geom_point(data = iq_90_ols_con[iq_90_ols_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "green4") +
  geom_errorbar(data = iq_90_ols_con[iq_90_ols_con$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "green4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.02, 0.0275) + 
  annotate(geom = "text", x = 65:95, y = -0.0165,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

temperature_deviation_iq_nb_con <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$median > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week.nb <- glm(constrained ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                              constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                              unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                              as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint_robust(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint_robust(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_iq[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_iq[1,], type = "response")) /  mean(iq_sig1$constrained)
    
    reg_temp_month.nb <- glm(constrained ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                               temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                               unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint_robust(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint_robust(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(constrained ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                                   as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint_robust(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint_robust(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(constrained ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                                    as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint_robust(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint_robust(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_nb_con <- temperature_deviation_iq_nb_con(iq_sig1, 0.95)
iq_90_nb_con <- temperature_deviation_iq_nb_con(iq_sig1, 0.90)

# Plot:

p_iq_nb_con <- ggplot() + 
  
  geom_line(data = iq_95_nb_con[iq_95_nb_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "green4") +
  geom_line(data = iq_95_nb_con[iq_95_nb_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_con[iq_95_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_con[iq_95_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_nb_con[iq_95_nb_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_con[iq_95_nb_con$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_con[iq_90_nb_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_con[iq_90_nb_con$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_con[iq_95_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_con[iq_95_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_con[iq_90_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_con[iq_90_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_con[iq_95_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_con[iq_95_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_con[iq_90_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_con[iq_90_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_con[iq_95_nb_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "green4") + 
  geom_errorbar(data = iq_95_nb_con[iq_95_nb_con$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "green4", linetype = "dashed") +
  geom_point(data = iq_90_nb_con[iq_90_nb_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "green4") +
  geom_errorbar(data = iq_90_nb_con[iq_90_nb_con$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "green4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.035, 0.05) + 
  annotate(geom = "text", x = 65:95, y = -0.03,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

# IED Attacks:
# OLS model results:
temperature_deviation_iq_ols_ied <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$median > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(log(ied+1) ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                            constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                            unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            week + governorate | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(ied+1) ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                             unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + governorate | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(ied+1) ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained |
                                 week + governorate | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(ied+1) ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained |
                                  month + governorate | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_ols_ied <- temperature_deviation_iq_ols_ied(iq_sig1, 0.95)
iq_90_ols_ied <- temperature_deviation_iq_ols_ied(iq_sig1, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
iq_95_ols_ied$change <- (sd(iq_sig1$median_dev, na.rm = T) * iq_95_ols_ied$estimate) / mean(log(iq_sig1$ied+1))

# Plot:
p_iq_ols_ied <- ggplot() + 
  
  geom_line(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "firebrick") +
  geom_line(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") + 
  geom_errorbar(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "firebrick", linetype = "dashed") +
  geom_point(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") +
  geom_errorbar(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "firebrick") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.02, 0.0275) + 
  annotate(geom = "text", x = 65:95, y = -0.0165,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

temperature_deviation_iq_nb_ied <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$median > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week.nb <- glm(ied ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                              constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                              unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                              as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint_robust(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint_robust(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_iq[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_iq[1,], type = "response")) /  mean(iq_sig1$ied)
    
    reg_temp_month.nb <- glm(ied ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                               temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                               unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint_robust(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint_robust(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(ied ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                                   as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint_robust(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint_robust(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(ied ~ median_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                                    as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "median_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$median, na.rm = T)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint_robust(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint_robust(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$median)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
  
}

iq_95_nb_ied <- temperature_deviation_iq_nb_ied(iq_sig1, 0.95)
iq_90_nb_ied <- temperature_deviation_iq_nb_ied(iq_sig1, 0.90)

# Plot:

p_iq_nb_ied <- ggplot() + 
  
  geom_line(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "firebrick") +
  geom_line(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") + 
  geom_errorbar(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "firebrick", linetype = "dashed") +
  geom_point(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") +
  geom_errorbar(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "firebrick") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.035, 0.05) + 
  annotate(geom = "text", x = 65:95, y = -0.03,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

iraq_combined <- ggarrange(p_iq_nb + xlab("") + ylab("Quasi-Poisson"), 
                           p_iq_nb_con + xlab("") + ylab(""), 
                           p_iq_nb_ied + xlab("") + ylab(""),  
                           p_iq_ols + xlab("") + ylab("OLS"),
                           p_iq_ols_con + xlab("") + ylab(""), 
                           p_iq_ols_ied + xlab("") + ylab(""), 
                           ncol = 3, nrow=2, labels  = "AUTO")

iraq_combined <- annotate_figure(iraq_combined,
                top = text_grob(expression(rho), size = 18, face = "bold"),
                bottom = text_grob(expression(gamma*" (°F)"), 
                                   face = "bold", size = 18),
                left = text_grob(expression(beta~Delta*"Temperature"), 
                                 size = 18, rot = 90, face = "bold"))

###
### Response results:
###

# Unconstrained Violence:
# QP model results:
iq_sig1$median2 <- iq_sig1$median^2
iq_sig1$median3 <- iq_sig1$median^3

IQ.w.u.l.qp <- bayesglm(unconstrained ~ median + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.u.q.qp <- bayesglm(unconstrained ~ median + median2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.u.c.qp <- bayesglm(unconstrained ~ median + median2 + median3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

IQ.m.u.l.qp <- bayesglm(unconstrained ~ median + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.u.q.qp <- bayesglm(unconstrained ~ median + median2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.u.c.qp <- bayesglm(unconstrained ~ median + median2 + median3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

# Constrained Violence:
IQ.w.c.l.qp <- bayesglm(constrained ~ median + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.c.q.qp <- bayesglm(constrained ~ median + median2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.c.c.qp <- bayesglm(constrained ~ median + median2 + median3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

IQ.m.c.l.qp <- bayesglm(constrained ~ median + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.c.q.qp <- bayesglm(constrained ~ median + median2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.c.c.qp <- bayesglm(constrained ~ median + median2 + median3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

# ied:
IQ.w.i.l.qp <- bayesglm(ied ~ median + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.i.q.qp <- bayesglm(ied ~ median + median2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.i.c.qp <- bayesglm(ied ~ median + median2 + median3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

IQ.m.i.l.qp <- bayesglm(ied ~ median + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.i.q.qp <- bayesglm(ied ~ median + median2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.i.c.qp <- bayesglm(ied ~ median + median2 + median3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

# Next, generate predicted counts:
# Again, construct a function for generating counts and creating output that can be plotted.
plot_function <- function(reg_output, X, dataset, type){
  set.seed(1)
  beta.sim <- mvrnorm(100, mu = coef(reg_output), Sigma = vcov(reg_output))
  paste <- model.matrix(reg_output)
  names1 <- colnames(paste)
  beta.sim <- as.data.frame(beta.sim)
  beta.sim <- beta.sim[, c(names1)]
  
  reps.df1 <- matrix(NA, length(na.omit(unique(round(dataset$median)))), 100)
  for(i in 1:length(na.omit(unique(round(dataset$median))))){
    if(X %in% "linear"){
      paste[, 2] <- i + (min(round(dataset$median), na.rm = T)-1)
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta)
      reps.df1[i, ] <- apply(yhat, 2, mean, na.rm=T)
    } else if(X %in% "quadratic"){
      
      paste[, 2] <- i + (min(round(dataset$median), na.rm = T)-1)
      paste[, 3] <- (i + (min(round(dataset$median), na.rm = T)-1))^2
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta)

      reps.df1[i, ] <- apply(yhat, 2, mean, na.rm=T)
    } else if(X %in% "cubic"){
      paste[, 2] <- i + (min(round(dataset$median), na.rm = T)-1)
      paste[, 3] <- (i + (min(round(dataset$median), na.rm = T)-1))^2
      paste[, 4] <- (i + (min(round(dataset$median), na.rm = T)-1))^3
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta)
      reps.df1[i, ] <- apply(yhat, 2, mean, na.rm=T)
    }
  }
  
  plot_dataset <- as.data.frame(matrix(NA, length(na.omit(unique(round(dataset$median)))),3))
  colnames(plot_dataset) <- c("means", "cil", "ciu")
  
  plot_dataset$means <- apply(reps.df1, 1 , mean)
  se <- apply(reps.df1, 1 , sd)
  plot_dataset$cil <- plot_dataset$means + qnorm(0.975)*se
  plot_dataset$ciu <- plot_dataset$means - qnorm(0.975)*se
  plot_dataset$temperature <- sort(na.omit(unique(round(dataset$median))))
  
  if(X %in% "linear"){
    plot_dataset$model <- "linear"
  } else if(X %in% "quadratic"){
    plot_dataset$model <- "quadratic"
  } else if(X %in% "cubic"){
    plot_dataset$model <- "cubic"
  }
  
  plot_dataset$max <- 0
  plot_dataset[which.max(plot_dataset$means), ]$max <- 1
  plot_dataset$attack_type <- type
  
  return(plot_dataset)
}

# Iraq:

p.IQ.w.u.l.qp <- plot_function(IQ.w.u.l.qp, "linear", iq_sig1, "unconstrained")
p.IQ.w.u.q.qp <- plot_function(IQ.w.u.q.qp, "quadratic", iq_sig1, "unconstrained")
p.IQ.w.u.c.qp <- plot_function(IQ.w.u.c.qp, "cubic", iq_sig1, "unconstrained")

p.IQ.m.u.l.qp <- plot_function(IQ.m.u.l.qp, "linear", iq_sig1, "unconstrained")
p.IQ.m.u.q.qp <- plot_function(IQ.m.u.q.qp, "quadratic", iq_sig1, "unconstrained")
p.IQ.m.u.c.qp <- plot_function(IQ.m.u.c.qp, "cubic", iq_sig1, "unconstrained")

p.IQ.w.c.l.qp <- plot_function(IQ.w.c.l.qp, "linear", iq_sig1, "constrained")
p.IQ.w.c.q.qp <- plot_function(IQ.w.c.q.qp, "quadratic", iq_sig1, "constrained")
p.IQ.w.c.c.qp <- plot_function(IQ.w.c.c.qp, "cubic", iq_sig1, "constrained")

p.IQ.m.c.l.qp <- plot_function(IQ.m.c.l.qp, "linear", iq_sig1, "constrained")
p.IQ.m.c.q.qp <- plot_function(IQ.m.c.q.qp, "quadratic", iq_sig1, "constrained")
p.IQ.m.c.c.qp <- plot_function(IQ.m.c.c.qp, "cubic", iq_sig1, "constrained")

p.IQ.w.i.l.qp <- plot_function(IQ.w.i.l.qp, "linear", iq_sig1, "ied")
p.IQ.w.i.q.qp <- plot_function(IQ.w.i.q.qp, "quadratic", iq_sig1, "ied")
p.IQ.w.i.c.qp <- plot_function(IQ.w.i.c.qp, "cubic", iq_sig1, "ied")

p.IQ.m.i.l.qp <- plot_function(IQ.m.i.l.qp, "linear", iq_sig1, "ied")
p.IQ.m.i.q.qp <- plot_function(IQ.m.i.q.qp, "quadratic", iq_sig1, "ied")
p.IQ.m.i.c.qp <- plot_function(IQ.m.i.c.qp, "cubic", iq_sig1, "ied")

# Plot:

# Iraq:
# Least Constrained:
plot.IQ.m.u.l.qp <- ggplot() + 
  geom_line(data = p.IQ.w.u.q.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.u.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  geom_line(data = p.IQ.w.u.c.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.u.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  
  geom_line(data = p.IQ.m.u.l.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) +
  geom_ribbon(data = p.IQ.m.u.l.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.IQ.m.u.q.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.IQ.m.u.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.IQ.m.u.c.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) +
  geom_ribbon(data = p.IQ.m.u.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(iq_sig1$temp)), max(round(iq_sig1$temp)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(2, 10) +
  ylab("") + 
  xlab("") +
  ggtitle("Least Constrained")

# Most Constrained:
plot.IQ.m.c.l.qp <- ggplot() + 
  geom_line(data = p.IQ.w.c.l.qp, aes(y = means, x = temperature)) +
  geom_ribbon(data = p.IQ.w.c.l.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "green4") +
  geom_line(data = p.IQ.w.c.q.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.c.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "green4") +
  geom_line(data = p.IQ.w.c.c.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.c.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "green4") +
  
  geom_line(data = p.IQ.m.c.l.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) +
  geom_ribbon(data = p.IQ.m.c.l.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.IQ.m.c.q.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.IQ.m.c.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.IQ.m.c.c.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) +
  geom_ribbon(data = p.IQ.m.c.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(iq_sig1$temp)), max(round(iq_sig1$temp)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(-0.5, 7.5) +
  ylab("") + 
  xlab("") +
  ggtitle("Most Constrained")

# IED:
plot.IQ.m.i.l.qp <- ggplot() + 
  geom_line(data = p.IQ.w.i.q.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.i.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "firebrick") +
  geom_line(data = p.IQ.w.i.c.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.i.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "firebrick") +
  
  geom_line(data = p.IQ.m.i.q.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.IQ.m.i.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.IQ.m.i.c.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.IQ.m.i.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(iq_sig1$temp)), max(round(iq_sig1$temp)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(5, 13) +
  ylab("") + 
  xlab("") +
  ggtitle("IED Attacks")

IQ.plots <- ggarrange(plot.IQ.m.u.l.qp + ylab(""), plot.IQ.m.c.l.qp, plot.IQ.m.i.l.qp, ncol = 3)
IQ.plots <- annotate_figure(IQ.plots,
                      bottom = text_grob(expression("Temperature (°F)"), 
                                         face = "bold", size = 18),
                      left = text_grob("Expected Count", 
                                       size = 18, rot = 90))
# PRODUCES FIGURE 16 IN PAPER:
ggarrange(IQ.plots, iraq_combined, ncol = 2)

###
### END.
###
